In [1]:
# importing libraries that I will be using.
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
import numpy as np
from sklearn.preprocessing import MinMaxScaler
import statsmodels.api as sm
import plotly.express as px
Border Crossing Dataset:¶
In [2]:
# reading in Border_Crossing_Entry_Data.csv and creating a dataframe.
border_crossing = pd.read_csv('Border_Crossing_Entry_Data.csv')
border_crossing
Out[2]:
Let's start getting to know the data! This will be able to tell us what we're working with, what we will and will not need and if there will be any cleaning necessary. I do a dataset at a time.
In [3]:
border_crossing.info()
# Date column will need to be changed to datetime. Latitude, Longitude, and Point can be dropped as they will not be needed (we
# have port codes as well as the port name, state and border entry area).
In [4]:
# Finding the number of duplicates.
border_crossing.duplicated().sum()
Out[4]:
In [5]:
# Finding the duplicates. Will need these later to drop them.
border_dupes = border_crossing.duplicated()
border_dupes2 = border_crossing[border_dupes]
border_dupes2
Out[5]:
As we can see above, these are not actual duplicates. We will still include these in the dataset.
In [6]:
# Finding the sum of all null values in the dataset.
border_crossing.isnull().sum()
Out[6]:
No need to worry about the nulls here as those columns will be dropped.
In [7]:
border_crossing['State'].unique()
Out[7]:
In [8]:
border_crossing['Port Name'].unique()
Out[8]:
In [9]:
border_crossing['Measure'].unique()
Out[9]:
In [10]:
border_crossing.nsmallest(8, 'Value')
Out[10]:
In [11]:
border_crossing.nlargest(8, 'Value')
Out[11]:
In [12]:
# Checking dates. These will need to match up with the drug_overdose dataset. Rows will need to br dropped accordingly.
border_max = border_crossing['Date'].max()
border_min = border_crossing['Date'].min()
print(border_max)
print(border_min)
Drug Overdose Dataset:¶
In [13]:
# reading in VSRR_Provisional_Drug_Overdose_Death_Counts.csv and creating a dataframe.
drug_overdose = pd.read_csv('VSRR_Provisional_Drug_Overdose_Death_Counts.csv')
drug_overdose
Out[13]:
In [14]:
drug_overdose.sample(8)
Out[14]:
In [15]:
drug_overdose.info()
In [16]:
# Checking for the years in this dataset. The border_crossing dataset goes back to 1997. If this doesn't go back as far, will
# isolate the years in the border_crossing set.
drug_overdose['Year'].unique()
Out[16]:
In [17]:
# Finding null values.
drug_overdose.isnull().sum()
Out[17]:
Data Cleaning:¶
border_crossing:
- Change the date column to datetime.
- Drop columns Latitude, Longitude, and Point.
- Drop all rows where 'Year' != 2015 and above. We are isolating what matches the drug_overdose dataset.
- Drop all rows where date is above 2023-08-01. This will match the drug_overdose set.
- Find ports with the highest values in each state. I only want to use ports in each state with the highest values.
- 'Melt' the 'Measure' values into one value. After trying multiple ways to use individual vehicles in the analysis, I've come to the conclusion that I won't be able to. I'm going to combine individual vehicles per month per state and port into one value.
- Change the data type for the 'Vehicle Total' from float64 to int64.
drug_overdose:
- Drop any 'indicator' rows that are not 'Number of Drug Overdose Deaths'.
- Combine 'Year' and 'Month' to create one datetime column.
- Figure out what YC is in the state column???
- Drop the 'State' Column. Rename 'State Name' to 'State'.
- Drop 'Percent Complete', 'Percent Pending Investigation', 'Footnote', 'Footnote symbol', and 'Predicted Value'.
- Deal with the nulls in 'Data Value'. I have a feeling this is because of certain rows that I won't be using anyway so once those are gone then I will reevaluate.
- Drop states that are not on the border_crossing datatset.
I will be starting with cleaning of the border_crossing dataset and moving on to the drug_overdose after that. This will ensure that things stay nice and tidy.¶
border_crossing cleaning:¶
In [18]:
# Change the date column to datetime.
border_crossing['Date'] = pd.to_datetime(border_crossing['Date'])
border_crossing.sample(5)
Out[18]:
In [19]:
# Making sure the 'Date' column was changed to datetime.
print(border_crossing.dtypes)
In [20]:
# Dropping Port Code, Latitude, Longitude, and Point columns as they are not needed. Keeping inplace=True because I just want to modify
# the original dataframe and not create a new one.
border_crossing.drop(columns=['Port Code', 'Latitude', 'Longitude', 'Point'], inplace=True)
border_crossing.sample(5)
Out[20]:
In [21]:
# Dropping all rows where year in 'Date' is below 2015.
border_crossing = border_crossing[border_crossing['Date'].dt.year >= 2015]
border_crossing.sample(5)
Out[21]:
In [22]:
# Checking to make sure there are no years below 2015 left.
checking_years = border_crossing['Date'].dt.year.unique()
print(checking_years)
In [23]:
# Checking dates. These will need to match with the drug_overdose dataset. Will drop rows accordingly.
border_min = border_crossing['Date'].min()
border_max = border_crossing['Date'].max()
print(border_min)
print(border_max)
In [24]:
# Filtering and dropping any rows with dates above 2023-08-01 as the drug_overdose dataset only goes to 2023-08-01.
border_dates = border_crossing['Date'] <= '2023-08-01'
border_filtered = border_crossing[border_dates]
border_filtered.sample(5)
Out[24]:
In [25]:
# Making sure the correct dates were dropped. Everything above 2023-08-01.
border_min2 = border_filtered['Date'].min()
border_max2 = border_filtered['Date'].max()
print(border_min2)
print(border_max2)
In [26]:
border_filtered['Measure'].unique()
Out[26]:
In [27]:
border_filtered['State'].unique()
Out[27]:
Finding the ports per state with the highest values.¶
These are the ports I will be using going forward. Doing multiple ports in each state and multiple types of vehicles was not working.
In [28]:
# Filtering and finding the rows by 'Port Name' and 'Value'. This will sort by the highest value by port alphabetically.
high_try = border_filtered.loc[border_filtered.groupby('Port Name')['Value'].idxmax()].sort_values(by='State')
# I only want to see these columns.
high_try_columns = high_try[['Port Name', 'State', 'Value']]
# I need to see all of the rows. This just tells pandas to display all the rows here.
pd.set_option('display.max_rows', None)
# Print it.
print(high_try_columns)
Ports with the highest values.¶
Ports that are not listed here will be filtered out and dropped.
- Alaska / Skagway
- Arizona / Nogales
- California / San Ysidro
- Idaho / Eastport
- Maine / Calais
- Michigan / Detroit
- Minnesota / International Falls
- Montana / Piegan
- New Mexico / Santa Teresa
- New York / Buffalo Niagara Falls
- North Dakota / Pembina
- Texas / El Paso
- Vermont / Highgate Springs
- Washington / Blaine
In [29]:
# Filtering border_filtered by the ports listed because it is faster to filter than drop over 100 different ports.
ports = ['Skagway', 'Nogales', 'San Ysidro', 'Eastport', 'Calais', 'Detroit', 'International Falls', 'Piegan', 'Santa Teresa',
'Buffalo Niagara Falls', 'Pembina', 'El Paso', 'Highgate Springs', 'Blaine']
border_final = border_filtered[border_filtered['Port Name'].isin(ports)]
# Checking to make sure the only ports are the ones listed above.
border_final['Port Name'].unique()
Out[29]:
Using pivot_table() to reshape the dataframe so that i can add all of the vehicle columns into one.¶
In [30]:
border_pivot = pd.pivot_table(border_final, values='Value', columns='Measure', index=['State', 'Date'], aggfunc="sum")
border_pivot.sample(5)
Out[30]:
Explaining replacing NaN values instead of imputing:¶
Originally when I used pivot_table() to reshape the dataframe to add all the vehicle columns, there ended up being many NaN values. I didn't understand why until I started locating the rows, specifically on 'State' and 'Date'. For some states, all vehicle types are recorded with associating values. However, there are some states that don't have certain vehicle types recorded at all.
I decided against imputing because I didn't want to add data with kNN that could influence analytics and models later. I simply decided to replace 'NaN' with pd.NA, this signifies that the data itself is actually missing after the pivot table adjustment.
Examples below.
In [31]:
# Example showing how Michigan has all vehicle types.
row_check = border_final.loc[(border_final['State'] == 'Michigan') & (border_final['Date'] == '2023-07-01')]
print(row_check)
In [32]:
# Example showing that California does not have all vehicle types.
row_check = border_final.loc[(border_final['State'] == 'California') & (border_final['Date'] == '2023-07-01')]
print(row_check)
In [33]:
border_pivot.replace(np.nan, pd.NA, inplace=True)
border_pivot.sample(5)
Out[33]:
In [34]:
# Creating a new column, 'Vehicle Total', by adding all of the vehicle type values into one.
border_pivot['Vehicle Total'] = border_pivot.sum(axis=1)
border_pivot.sample(5)
Out[34]:
In [35]:
border_pivot.drop(columns=['Bus Passengers', 'Buses', 'Pedestrians', 'Personal Vehicle Passengers', 'Personal Vehicles',
'Rail Containers Empty', 'Rail Containers Loaded', 'Train Passengers', 'Truck Containers Empty',
'Truck Containers Loaded', 'Trucks', 'Trains'], inplace=True)
border_pivot.sample(5)
Out[35]:
In [36]:
border_pivot.info()
In [37]:
# Changing the data type to int64.
border_pivot['Vehicle Total'] = border_pivot['Vehicle Total'].astype('int64')
border_pivot.info()
drug_overdose cleaning:¶
In [38]:
# Dropping the 'State' column.
drug_overdose.drop(columns=['State'], inplace=True)
drug_overdose.sample(5)
Out[38]:
In [39]:
# Renaming 'State Name' to 'State'.
drug_overdose = drug_overdose.rename(columns={'State Name':'State'})
drug_overdose.sample(5)
Out[39]:
Combining the 'Year' and 'Month' column into one column called 'Date'. This is going to require bringing together strings with a hyphen and then converting to datetime. Before that can happen, the 'Year' column has to be converted to an object.¶
In [40]:
# Converting 'Year' column to an object.
drug_overdose['Year'] = drug_overdose['Year'].astype(str)
print(drug_overdose.dtypes)
In [41]:
drug_overdose.info()
In [42]:
# Combining the 'Year' and 'Month' with a hyphen.
drug_overdose['Date'] = drug_overdose['Year'].astype(str) + '-' + drug_overdose['Month'].astype(str)
drug_overdose.sample(5)
Out[42]:
In [43]:
# Converting 'Date' over to datetime.
drug_overdose['Date'] = pd.to_datetime(drug_overdose['Date'])
drug_overdose.sample(5)
Out[43]:
In [44]:
# Dropping 'Period', 'Percent Complete', 'Percent Pending Investigation', 'Footnote', 'Footnote symbol', 'Predicted Value',
# 'Year' and 'Month' as I won't be needing them anymore.
drug_overdose.drop(columns=['Period', 'Percent Complete', 'Percent Pending Investigation', 'Footnote', 'Footnote Symbol', 'Predicted Value', 'Year', 'Month'], inplace= True)
drug_overdose.sample(5)
Out[44]:
In [45]:
# Finding unique values within 'Indicator' that I will be dropping except, 'Number of Drug Overdose Deaths'.
drug_overdose['Indicator'].unique()
Out[45]:
In [46]:
# Filtering and dropping rows that do not contain, 'Number of Drug Overdose Deaths'.
drug_filter = drug_overdose['Indicator'] == 'Number of Drug Overdose Deaths'
drug_overdose = drug_overdose[drug_filter]
drug_overdose.sample(5)
Out[46]:
In [47]:
# Making sure that the null values have been dealt with. What I assumed was correct, they were in catergories I didnt' need.
drug_overdose.isnull().sum()
Out[47]:
In [48]:
# Checking dates. These will need to match with the border_crossing dataset. Will drop rows accordingly.
drug_max = drug_overdose['Date'].max()
drug_min = drug_overdose['Date'].min()
print(drug_max)
print(drug_min)
Dropping the states that are not border states is next. Dropping rows based on specific column ('State') values. They include: Maine, Texas, Washington, North Dakota, Arizona, Vermont, Montana, Alaska, New York, Michigan, Minnesota, New Mexico, California, and Idaho.¶
In [49]:
drop_states = [ 'Alabama', 'Arkansas', 'Colorado', 'Connecticut', 'Delaware', 'District of Columbia', 'Florida', 'Georgia',
'Hawaii', 'Illinois', 'Indiana', 'Iowa', 'Kansas', 'Kentucky', 'Louisiana',
'Maryland', 'Massachusetts', 'Mississippi', 'Missouri', 'Nebraska', 'Nevada',
'New Hampshire', 'New Jersey', 'New York City', 'North Carolina', 'Ohio', 'Oklahoma', 'Oregon',
'Pennsylvania', 'Rhode Island', 'South Carolina', 'South Dakota', 'Tennessee',
'Utah', 'Virginia', 'West Virginia', 'Wisconsin','Wyoming', 'United States']
drug_overdose = drug_overdose[drug_overdose.State.isin(drop_states) == False]
drug_overdose['State'].unique()
Out[49]:
In [50]:
# Making sure that the states listed above match the border_filtered dataset.
border_sorted = sorted(border_filtered['State'].unique())
border_sorted
Out[50]:
In [51]:
# Will no longer be needing 'Indicator' column as, 'Data Value' column is synonomus with it.
drug_overdose.drop(columns=['Indicator'], inplace= True)
drug_overdose.sample(5)
Out[51]:
In [52]:
# Renaming'Data Value' with 'Drug Overdose Deaths'. This helps keeps things tidy when merging the datasets.
drug_overdose = drug_overdose.rename(columns={'Data Value': 'Drug Overdose Deaths'})
drug_overdose.sample(5)
Out[52]:
In [53]:
# Checking to make sure data types are correct. I will be changing 'Drug Overdose Deaths' column to int64.
drug_overdose.info()
In [54]:
# I keep running into an error when trying to change the 'Drug Overdose Deaths' column into int64. The error has
# to do with there being commas in the values. Now I'm curious how many there actually are.
comma_count = drug_overdose['Drug Overdose Deaths'].str.contains(',').sum()
print(comma_count)
In [55]:
# Now that I know there a lot of commas, I'm going to get rid of them. They are not needed and I need to convert
# the 'Drug Overdose Deaths' column.
drug_overdose['Drug Overdose Deaths'] = drug_overdose['Drug Overdose Deaths'].str.replace(',', '')
drug_overdose['Drug Overdose Deaths'] = drug_overdose['Drug Overdose Deaths'].astype('int64')
drug_overdose.info()
Combining the datasets together:¶
In [56]:
final_border = pd.merge(border_pivot, drug_overdose, on=['State', 'Date'])
final_border.sample(10)
Out[56]:
In [57]:
# Making sure that merge() worked and there are no NaN values. Success!
final_border.isna().sum()
Out[57]:
In [58]:
final_border.info()
In [59]:
final_check = final_border.loc[(final_border['State'] == 'Texas') & (final_border['Date'] == '2018-07-01')]
print(final_check)
1. Visualizing descriptive statistics¶
Overview of the distribution of drug overdose deaths.¶
For the descriptive statistics part of this project I wanted to really understand the spread of what drug overdoses and vehicle frequency look like in states that border the U.S. Mexico and U.S. Canada border. You will see that they are broken up by drug overdoses and vehicle frequency and then what they look like together. I started out using describe() on the 'Drug Overdose Deaths' column, which gives a basic overview of elementary statistics. The average amount of drug overdose deaths per month is around ~1501 and the max was ~12000. Below are more visualizations that will help to understand what drug overdose deaths in border states look like.
In [60]:
final_border['Drug Overdose Deaths'].describe()
Out[60]:
The graph below is a Kernel Density Plot via Seaborn. Originally I had used pandas own plot() and added multiple components such as 'kind', 'xlim', and 'xticks'. While these two graphs were very similar, I found that using seaborn kdeplot() worked much better and elimated the need for extra elements that plot() requires. I can see that drug overdose deaths were congrugated below 1000 per month and steadily decrease after 1000. This is supported by the statistics above. There are also outliers around the 12000 mark.¶
Kernel Density Plot can be thought of as a much smoother version of a histogram. Instead of needing to put data into bins, it uses a smoothing function below the hood to visualize density distrubution.¶
In [61]:
DrugOD_column = final_border['Drug Overdose Deaths']
plt.figure(figsize=(17,8))
sns.kdeplot(DrugOD_column)
Out[61]:
The graph below is a visualization of drug overdose deaths per year, 2015-mid 2023. I used lineplot() via seaborn to show the mean and confidence interval over time. The dark line on the graph represents the mean drug overdose deaths per year, the shaded region behind the line is a 95% confidence interval. The confidence interval represents the range of possible values around the true mean. We can see a uptick in drug overdose deaths during and after 2020.¶
The issue with this graph is that is doesn't respresent the entirety of drug overdose deaths over time. In the Kernel Density graph above, we can see that even though the deaths above 2500 per month were small, they were still there. To understand that, I used a box graph next.¶
In [62]:
plt.figure(figsize=(17,7))
sns.lineplot(x='Date', y='Drug Overdose Deaths', data=final_border)
Out[62]:
I utilized the boxplot() via seaborn to visualize and understand the outliers in drug overdose deaths per year. The line in the middle of the boxplot is the median. The box show how spread out around the median the data is per year.¶
In [63]:
plt.figure(figsize=(17,9))
final_border['Year'] = final_border['Date'].dt.year
sns.boxplot(x='Year', y='Drug Overdose Deaths', width=0.5, data=final_border)
Out[63]:
In [64]:
plt.figure(figsize=(17,8))
sns.barplot(x='State', y='Drug Overdose Deaths', data=final_border)
Out[64]:
Overview of entering vehicle frequency at the U.S. Mexico and U.S. Canada border.¶
In [65]:
final_border['Vehicle Total'].describe()
Out[65]:
The below graph is a Kernel Density plot. I chose this because it's smoother than a histogram, instead of having to put values into bins, as well as eliminating the need for extra elements to be added. The x-axis is by millions, we can see that the majority of the vehicles that pass through monthly are around 800,000, this would be on par with the statistics above. It tapers off around 1.5 million and experiences a small bump around 4 million per month. If I were to overlay this graph and the drug overdoses per month together, they would be very similar.¶
In [66]:
vehicle_freq = final_border['Vehicle Total']
plt.figure(figsize=(17,8))
sns.kdeplot(vehicle_freq)
Out[66]:
This bar plot below shows total vehicle frequency by state. We can see that the majority of vehicle frequency is passing through the U.S. Mexico border and less so on the U.S. Canada border, the largest vehicle frequency being in California and Texas.¶
In [67]:
plt.figure(figsize=(17,8))
plt.xticks(rotation=45)
sns.barplot(x='State', y='Vehicle Total', data=final_border)
Out[67]:
The line plot below shows us that vehicle frequency was consitent from 2015-2019. During the summer months vehicle frequency would peak and during the winter months, it would dip back down. The consitency lasted until early 2020 when the Covid-19 lockdowns happened. During this time the borders to both Mexico and Canda were locked down. Vehicle freqeuncy has slowly picked back up again, around mid 2022 is when vehicle frequency returned to pre Covid-19 levels.¶
In [68]:
plt.figure(figsize=(17,7))
sns.lineplot(x='Date', y='Vehicle Total', data=final_border)
Out[68]:
Boxenplots are are great way to visualize outliers in state vehicle frequency.¶
In [69]:
plt.figure(figsize=(17,9))
sns.boxenplot(data=final_border, x="Vehicle Total", y="State")
Out[69]:
Overview of drug overdose deaths and entering vehicle frequency at the U.S. Mexico and U.S. Canada border together.¶
Scaling is going to be required for the two columns that I want to visualize. This is because the average value from 'Vehicle Total' is 800,000 while the average value from 'Drug Overdose Deaths' column is 1599. MinMaxScaler scales data that is different so that it can be visualized and modeled correctly. The scaled dataframe will also be used for correlation and regression later.¶
In [70]:
# Defining the scaler.
scaler = MinMaxScaler()
# Fitting the scaler and transforming two columns.
final_border[['Vehicle Total', 'Drug Overdose Deaths']] = scaler.fit_transform(final_border[['Vehicle Total',
'Drug Overdose Deaths']])
final_border.sample(5)
Out[70]:
After the data has been scaled for both Vehicle Total and Drug Overdose Death columns, I wanted to visualize them. Below is scaled kde plot both columns, similar to the individual plots above.¶
In [71]:
plt.figure(figsize=(17,7))
scaled_DrugOD_column = final_border['Drug Overdose Deaths']
scaled_vehicle_total = final_border['Vehicle Total']
plt.figure(figsize=(17,8))
sns.kdeplot(scaled_DrugOD_column, fill='yes')
sns.kdeplot(scaled_vehicle_total, fill='yes')
Out[71]:
Visualizing both Drug Overdose Deaths and Vehicle Freqeuncy over time, from 2015-mid 2023¶
In [72]:
plt.figure(figsize=(17,7))
# Using line plot and getting rid of the confidence Interval.
sns.set_style('darkgrid')
sns.lineplot(x='Date', y='Drug Overdose Deaths', data=final_border, errorbar=None, color='green')
sns.lineplot(x='Date', y='Vehicle Total', data=final_border, errorbar=None, color='orange')
Out[72]:
Visualizing a scatterplot by state for both drug overdose deaths and vehicle freqeuncy.¶
In [73]:
plt.figure(figsize=(17,9))
sns.scatterplot(data=final_border, x="Drug Overdose Deaths", y="Vehicle Total", hue="State")
Out[73]:
Left bottom corner zoomed in version of the plot above.¶
In [74]:
plt.figure(figsize=(17,9))
sns.scatterplot(data=final_border, x="Drug Overdose Deaths", y="Vehicle Total", hue="State")
plt.xlim(0.0,0.4)
plt.ylim(0.0,0.4)
Out[74]:
2. Visualizing Correlation¶
Correlation is the measurement of the strength of the relationship between two variables. Below I found correlation between drug overdose deaths and vehicle frequency, on a scale of -1 (neg relationship) and 1 (positive relationship), the correlation here falls at 0.8! There is a positive relationship, when vehicle total begins to rise, drug overdose deaths do as well.
In [75]:
correlation = final_border['Vehicle Total'].corr(final_border['Drug Overdose Deaths'])
print(correlation)
In [76]:
sns.lmplot(data=final_border, x="Drug Overdose Deaths", y="Vehicle Total", height=7)
Out[76]:
Hypothesis Test¶
Null Hypothesis:
$H_{0}$: The less vehicles that pass through the most frequented ports in each state, the less occurrence of overdose deaths in U.S. Mexico and U.S. Canda bordering states there will be.
Alternative Hypothesis:
$H_{A}$: The more vehicles that pass through the most frequented ports in each state, the larger the occurrence of overdose deaths in U.S. Mexico and U.S. Canda bordering states there will be.
Planned Statistical Test:
Linear Regression
Metrics grenerated from planned statistical test:
T-statistic, P-Value
Alpha Value (threshold):
$a$: 0.05%
if p-value is < 0.05 we reject the null hypothesis.
if p-value is > 0.05 we fail to reject the null hypothesis.
In [77]:
# Defining the independent and dependent variables. Using only the only two columns that are needed here.
x = final_border['Vehicle Total']
y = final_border['Drug Overdose Deaths']
# Adding a constant.
x_const = sm.add_constant(x)
# fit the model
regression_model = sm.OLS(y, x_const).fit()
# printing summary statistics
print(regression_model.summary())
P-value is 0.05% we reject the null hypothesis. This means that there is a strong relationship between vehicle frequency and drugoverdose deaths. It is unlikely due to chance that there is a relationship. The alternative hypothesis is true.¶
In [78]:
# Trying plotly for the first time. I should have been using this all along.
fig = px.scatter(final_border, x="Drug Overdose Deaths", y="Vehicle Total", trendline="ols", trendline_color_override="orange")
fig.show()
In [ ]: